1.1 LaTeX

R Markdown est encore plus puissant, on peut utilisier \(\LaTeX\), soit inline comme \(y_i=\beta_0 + x_i\beta_1\) ou comme equation:

\[ \hat{\beta}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y} \] vous pouvez même inclure des fichiers .tex entiers.

1.2 Figures

Images (ce n’est pas un plot) peuvent également être facilement incorporés avec ![](la_meilleure_image.png)

1.3 Knitting

Pour les documents de base, le yaml est ce qui est important:

rmdformats::robobook c’est seulement pour rendre le ficher légèrement plus joli.

1.4 Pourquoi?

• Projets pour votre thèse de baccalauréat (Peut être entièrement reproductible - un click)
• après votre baccalauréat
  • Industrie-1: Vous ferez beaucoup de rapports, souvent répétitifs, ce qui n’est pas très fun.. mais vous pouvez automatiser la plupart avec markdown
  • Industrie-2: Particulierement en tant que junior, vous aurez besoin de faire des recherches, Rmarkdown est parfait pour les présenter
  • Université: Si vous faites un master - Vous devrez probablement coder et remettre votre code
  • Université: De plus en plus les “journals” exigent le code pour les publications
• Vous pouvez également utiliser python/bash/stan à partir des chunks de rmarkdown / à mon avis c’est mieux que jupyter notebook

2.1 SQL-Like

SQL (Structured Query Language) - votre meilleur ami si vous travaillez avec des données dans une entreprise, mais également utile à plus petite échelle (aka. pour vos projets).

L’exercice est basé sur l’article : “Where is the Land of Opportunity ? The Geography of Intergenerational Mobility in the United States” (Raj Chetty, Nathaniel Hendren, Patrick Kline et Emmanuel Saez ; Quarterly Journal of Economics 129(4) : 1553-1623, 2014).Il fait partie d’un plus grand projet disponible à l’adresse suivante : https://opportunityinsights.org/.

Nous commençons par charger quelques données échantillons. Notez que “kable” dans la dernière ligne rend seulement les graphiques jolis.

library(haven) # si le format est .dta
library(readxl) # se le format est .xls / .xlsx

mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")

mobility_data %>% 
  head() %>% 
  kable()

Nous voyons qu’il y a des données manquantes et que les variables explicatives sont normalement numériques. Nous nous en occuperons plus tard. Pour une explication des variables, considérez le fichier codebook.xls.

select

En général, un jeu de données comporte de nombreuses colonnes. Mais parfois, nous n’avons besoin de considérer que certaines d’entre elles. La logique “select” vous permet de choisir un sous-ensemble des colonnes de données.

Il existe de nombreuses façons de le faire. Ci-dessous la démonstration pour sélectionner une seule colonne “e_rank_b” et une collection de colonnes c(“e_rank_b”, “taxrate”, “crime_violent”). C’est quoi la difference?

# Une seule colonne
facon_1 <- mobility_data$e_rank_b
facon_2 <- mobility_data[,'e_rank_b']
facon_3 <- mobility_data %>% select(e_rank_b)

# Plusieurs
facon_1_mul <- mobility_data[, c("e_rank_b", "taxrate", "crime_violent")]
facon_2_mul <- mobility_data %>% select(e_rank_b, taxrate, crime_violent)


facon_1 %>% head() 

## [1] 43.58841 39.64703 46.93892       NA 43.22648 53.63568

facon_1_mul %>% head()

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     39.6  0.0257        0.0316
## 3     46.9  0.0117        0.0512
## 4     NA    0.0217        0.0371
## 5     43.2  0.0452        0.0500
## 6     53.6  0.0267        0.0322

facon_2 %>% head() 

## # A tibble: 6 × 1
##   e_rank_b
##      <dbl>
## 1     43.6
## 2     39.6
## 3     46.9
## 4     NA  
## 5     43.2
## 6     53.6

facon_2_mul %>% head() 

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     39.6  0.0257        0.0316
## 3     46.9  0.0117        0.0512
## 4     NA    0.0217        0.0371
## 5     43.2  0.0452        0.0500
## 6     53.6  0.0267        0.0322

if-condition (filter)

Les commandes de “select” agissent sur les colonnes, les équivalents pour les lignes sont les “if-conditions” dans R, on appel ca les conditions “filter”. En reprenant l’exemple ci-dessus. Nous ne voulons que les observations où le e_rank_b est supérieur à 40. Quelle différence peut-on constater entre les deux méthodes ?

facon_1 <- mobility_data[mobility_data$e_rank_b > 40, # En base-R on travaille avec 
                                                      # un vecteur logique
                         c("e_rank_b", "taxrate", "crime_violent")]

facon_2 <- mobility_data %>% 
  select(e_rank_b, taxrate, crime_violent) %>% 
  filter(e_rank_b > 40)

facon_1 %>% head()

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     46.9  0.0117        0.0512
## 3     NA   NA            NA     
## 4     43.2  0.0452        0.0500
## 5     53.6  0.0267        0.0322
## 6     45.6  0.0412        0.0267

facon_2 %>% head()

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     46.9  0.0117        0.0512
## 3     43.2  0.0452        0.0500
## 4     53.6  0.0267        0.0322
## 5     45.6  0.0412        0.0267
## 6     48.8  0.0222        0.0286

C’est important à garder dans l’esprit les differences entre tidyverse et base R. En practique, on jete les observations manquantes souvent, mais il faut toujours faire le “sanity check”!

groupby

Parfois, nous sommes intéressés par des données agrégées. Avec la commande group_by, nous pouvons facilement créer des statistiques sommaires pour des sous-groupes de données (et on peut les ajouter à la base originale).

Il existe des façons de faire exactement cela dans base R, mais ici nous nous concentrons sur la logique de dplyr, qui est plus intuitive.

# Calculer la moyenne de "e_rank_b" pour chaque ville

group_stat <- mobility_data %>% # Choisir jeu des données
  select(czname, e_rank_b) %>% # Choisir les colonnes
  group_by(czname) %>% # Specifier les groups
  summarise(
    mean_per_city = mean(e_rank_b, na.rm=T) # Specifier la statistique d'aggreagtion
  )

group_stat %>% head()

## # A tibble: 6 × 2
##   czname     mean_per_city
##   <chr>              <dbl>
## 1 Aberdeen            45.6
## 2 Aiken               38.5
## 3 Ainsworth           48.6
## 4 Albany              52.9
## 5 Alexandria          46.6
## 6 Allentown           61.3

(mutate)

Souvent, nous devons transformer les variables pour obtenir des informations significatives de notre analyse. Cela peut être fait facilement avec la logique “mutate”. L’utilisation de condtions if-else pour créer de nouvelles variables est extrêmement utile. Par exemple, nous pouvons diviser la mobilité sociale (e_rank_b) en deux catégories : faible et élevée.

data_tmp <- mobility_data
data_tmp <- data_tmp[!is.na(data_tmp$e_rank_b), ]
data_tmp$binary_mobility <- ifelse(data_tmp$e_rank_b > 45, 'High', 'Low')

ex_mutate_2 <- mobility_data %>% 
  filter(!is.na(e_rank_b)) %>% # PQ? 
  mutate(binary_mobility = if_else(e_rank_b > 45, 'High', 'Low')) %>% 
  select(binary_mobility, e_rank_b)


data_tmp[,c('binary_mobility', 'e_rank_b')] %>% head()

## # A tibble: 6 × 2
##   binary_mobility e_rank_b
##   <chr>              <dbl>
## 1 Low                 43.6
## 2 Low                 39.6
## 3 High                46.9
## 4 Low                 43.2
## 5 High                53.6
## 6 Low                 38.1

ex_mutate_2 %>% head()

## # A tibble: 6 × 2
##   binary_mobility e_rank_b
##   <chr>              <dbl>
## 1 Low                 43.6
## 2 Low                 39.6
## 3 High                46.9
## 4 Low                 43.2
## 5 High                53.6
## 6 Low                 38.1

2.2 Exercises

1. Transformer la variable “e_rank_b” en “inférieure à 50” et “supérieure à 50”. Calculez la moyenne de “crime_violent” pour chaque groupe. Quelles sont vos conclusions?
2. Transformer la variable “frac_middleclass” en “inférieure à 0.5” et “supérieure à 0.5”. Pareil pour la variable “eitc_exposure”. Calculer la moyenne du taux d’abandon scolaire au lycée “dropout_r” pour chaque combinaison des deux nouvelles variables. Conseil : utilisez group_by
3. Proposer une façon de traiter les valeurs “NA”. Pourquoi pensez-vous que c’est une bonne idée ?

3.1 Plots

Les plots/figures peuvent facilement être intégrées dans markdown. Tout comme une sortie de tableau, on peut les structurer légèrement. Par exemple, nous pouvons ajouter une légende ou ajouter plusieurs graphiques côte à côte.

Par exemple, commençons par un histogramme:

# Histogram simple
hist(mobility_data$e_rank_b)

?hist

?hist

plot(density(data_tmp$e_rank_b)) # Pourquoi changer les données?
plot(density(data_tmp$e_rank_b), col='red')

Il existe aussi un package sympa de l’écosystème tidyverse (ggplot2)

plot_1 <- mobility_data %>% 
  ggplot() + 
  geom_histogram(aes(x=e_rank_b))

plot_2 <- mobility_data %>% 
  ggplot() + 
  geom_histogram(aes(x=crime_violent))

grid.arrange(plot_1, plot_2, ncol=2)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 19 rows containing non-finite values (stat_bin).

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 17 rows containing non-finite values (stat_bin).

Nous pouvons également construire des figures avec plusieurs informations. Par exemple, un histogramme avec une ligne de densité. (Attention toutefois aux axes!)

# Base R
hist(data_tmp$e_rank_b, probability = T)
lines(density(data_tmp$e_rank_b), col='red')

# GGplot
data_tmp %>% 
  ggplot() + 
  geom_histogram(aes(x=e_rank_b, y=..density..)) + 
  geom_density(aes(x=e_rank_b), color='red')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.2 Exercises

Conseil : si vous n’arrivez pas à trouver la bonne commande pour une figure, les solutions, utilisez google!

1. Montrer la relation entre le taux d’obtention d’un diplôme universitaire et la mobilité socialeà l’aide d’un diagramme de dispersion (scatterplot) et d’une ligne de régression. Conseil: vous pouvez trouver une solution si vous googlez “r plot with regression line”

2. Existe-t-il une relation entre les dépenses des gouvernements locaux par habitant (subcty_total_expenditure_pc) et le taux d’imposition (taxrate)? Pouvez-vous utiliser les données directement ? Conseil : pour mettre les données à l’échelle, vous pouvez toujours essayer d’utiliser log, sqrt, etc.

3. Montrez graphiquement qu’il existe une ville où la criminalité violente est exceptionnellement élevée (‘crime_violent’). Que pouvez-vous faire avec de telles données (outlier ou pas)?

4.1 Tables

Les tableaux peuvent facilement être incorporés dans R. Vous pouvez soit les imprimer directement, soit utiliser kable (du paquet knitr - plus jolie). Bien sûr, nous pouvons également combiner ce que nous avons vu dans les sections précédentes (en particulier la manipulation des données)

# Impression de base
mobility_data %>% head()

## # A tibble: 6 × 39
##   index    cz czname         e_rank_b cs_race_bla cs_race_theil_20… cs00_seg_inc
##   <dbl> <dbl> <chr>             <dbl>       <dbl>             <dbl>        <dbl>
## 1     0   100 Johnson City       43.6    0.0694             0.0577        0.0409
## 2     1   200 Morristown         39.6    0.445              0.0750       -0.0259
## 3     2   301 Middlesborough     46.9    0.0151             0.110         0.0954
## 4     3   302 Knoxville          NA     -0.000411          -0.00446       0.0105
## 5     4   401 Winston-Salem      43.2   -0.0193             0.263         0.0289
## 6     5   402 Martinsville       53.6    0.00563            0.00699       0.0221
## # … with 32 more variables: cs00_seg_inc_pov25 <dbl>, cs00_seg_inc_aff75 <dbl>,
## #   frac_traveltime_lt15 <dbl>, hhinc00 <dbl>, gini <dbl>,
## #   inc_share_1perc <dbl>, gini99 <dbl>, frac_middleclass <dbl>, taxrate <dbl>,
## #   subcty_total_expenditure_pc <dbl>, tax_st_diff_top20 <dbl>,
## #   eitc_exposure <dbl>, ccd_exp_tot <dbl>, ccd_pup_tch_ratio <dbl>,
## #   score_r <dbl>, dropout_r <dbl>, num_inst_pc <dbl>, tuition <dbl>,
## #   gradrate_r <dbl>, cs_labforce <dbl>, cs_elf_ind_man <dbl>, …

# Avec Kable
mobility_data %>% head() %>% kable() # c'est comme kable(head(mobilitiy_data))

ici nous pouvons manipuler légèrement nos données

# Example avec manipulation des données 
# Décrire ce qui se passe ici
mobility_data %>% 
  filter(!is.na(e_rank_b)) %>% 
  mutate(decile = ntile(e_rank_b, 10)) %>% 
  select(decile, e_rank_b, mig_outflow) %>% 
  group_by(decile) %>% 
  summarise(average_outflow = mean(mig_outflow)) %>% 
  kable()

4.2 Distributions en R

L’un des exercices les plus importants lors de la modélisation statistique consiste à vérifier les hypothèses du modèle. La logique est toujours la même.

q-distribution (quantiles, ex: qnorm, qpois, qgamma etc.) Quel est le 99ème quantile d’une gamma avec shape=2 et rate=1?

qgamma(0.99, shape=2,rate=1)

## [1] 6.638352

d-distribution (densité, ex: dnorm, dpois) Tracez visuellement le 99ème quantile d’une gamma avec shape=2 et rate=1?

gamma_dist <- dgamma(seq(0,8,0.1), shape=2,rate=1)

plot(x=seq(0,8,0.1), 
     y=gamma_dist,
     type='l')
abline(v=qgamma(0.99, shape=2,rate=1), 
       col='red', lty=2)

p-distribution (densité cumulative, ex: pnorm, ppois)

Supposons que \(X \sim Gamma(2,1)\). Quel pourcentage de la masse totale se trouve à droite de 5 ?

1 - pgamma(5, shape=2, rate=1)

## [1] 0.04042768

r-distribution (échantillon, ex: rnorm, rpois)

Simulez 100 observations à partir d’un gamma avec shape=2, rate=1 et comparez-les à la distribution théorique.

set.seed(42)
sims_gamma <- rgamma(100, shape=2, rate=1)

plot(density(sims_gamma), ylim=c(0,0.38), 
     xlim=c(-1, 7), 
     main='Comparaison empirique vs. theoretique')
lines(x=seq(0,8,0.1), 
      y=gamma_dist, col='red')

# A la main
# Sans transformation
inc <- (as.numeric(mobility_data$hhinc00))
inc <- (inc - mean(inc)) / sd(inc)

quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(inc, seq(0,1,0.001))

plot(quantiles_norm,
     quantiles_empirical, 
     main='Normal Q-Q Plot', 
     ylab='Sample Quantiles',  
     xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")

# Avec transformation
# A la main
# Sans transformation
log_inc <- log(as.numeric(mobility_data$hhinc00))
log_inc <- (log_inc - mean(log_inc)) / sd(log_inc)

quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(log_inc, seq(0,1,0.001))

plot(quantiles_norm,
     quantiles_empirical, 
     main='Normal Q-Q Plot (Log transformed)', 
     ylab='Sample Quantiles',  
     xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")

qqnorm(log(as.numeric(mobility_data$hhinc00)))
qqline(log(as.numeric(mobility_data$hhinc00)),
       col = 2,lwd=2,lty=2)

4.3 Exercises

1. Supposons que \(X \sim \mathcal{N}(0,10)\). Echantillonnez 10 observations et calculez la moyenne, que nous appelons Y. Répétez l’opération ci-dessus 1000 fois et calculez la moyenne et l’écart-type de Y. A quoi vous attendez-vous et pourquoi, qu’observez-vous ?

2. Créez un tableau ou qui montre le nombre d’observations manquantes pour chaque variable (colonne dans l’ensemble de données).

3. Analysez les variables tuition et num_inst_pc et discutez de ce que vous voyez. Quelles sont vos conclusions ? S’agit-il d’un problème?

4. Créez un tableau montrant le coefficient de corrélation entre la mobilité sociale e_rank_b et toutes les autres variables. Évaluez vos résultats.

5. Séparez la variable en “faible”, “moyen” et “élevé”. Justifiez votre choix. Comparez chaque niveau avec les variables gini et cs_born_foreign. Interprétez et créez un tableau.

5.1 Exercise

Notre objectif est de mieux comprendre la mobilité intragénérationnelle. On dit que la mobilité intergénérationnelle est élevée si le revenu des parents (au cours de leur vie) n’est pas très lié au revenu de leurs enfants. Pour les exercices, vous pouvez utiliser n’importe quel “package” que vous jugez nécessaire mais vous devez justifier vos résultats.

Pour cela, résolvez les exercices suivants :

1. Chargez les données, regardez les types de données, est-ce que cela a du sens ?

2. La variable d’intérêt est e_rank_b, renommez cette colonne. Que se passe-t-il si vous utilisez un espace ” ” dans le nom ? Est-ce une bonne idée ?

3. Visualiser la distribution de la mobilité intergénérationnelle. Donnez à vos graphiques un titre approprié et des étiquettes d’axe correctes.

4. Examinez la relation entre la mobilité intergénérationnelle et la taille du secteur manufacturier. Quelle est votre conclusion ?

5. Découpez la variable de la taille du secteur manufacturier en trois composantes (low-mid-high). Calculez les moyennes par groupe pour trois autres variables que vous jugez intéressantes.

6. Ajustez une régression qui décrit l’équation suivante \(\mathbb{E}[y| \text{manufacturing}]=\hat{\beta}_0 + \hat{\beta}_1*\text{manuf_1} + \hat{\beta}_2*\text{manuf_2}\). Pourquoi exclure la dernière composante?

L’existence d’une classe moyenne importante est souvent considérée comme une bonne chose. Voir par exemple cet article, analysons cela plus en détail.

7. Réalisez un nuage de points (“scatterplot”) avec la mobilité intergénérationnelle sur l’axe des y et la fraction de la classe moyenne sur l’axe des x. Quelles sont vos conclusions ?

8. Exécutez un modèle linéaire qui explique la mobilité intergénérationnelle à l’aide de la fraction de la classe moyenne.

9. Nous voulons rendre le modèle légèrement plus flexible. Effectuez 10 régressions différentes qui correspondent à des polynômes de plus en plus complexes de la fraction de la classe moyenne.

10. Défi : Visualisez les différents résultats. Qu’arrive-t-il à la valeur \(R^2\)? Est-ce une bonne chose ?

5.2 Charger et preparer les données

1. Chargez les données, regardez les types de données, est-ce que cela a du sens ?

# Charger l'essentiel
library(tidyverse)
library(knitr)
library(gridExtra)
library(haven) 
library(readxl)

mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")

mobility_data %>% 
  head() %>% 
  kable()

str(mobility_data)

## tibble [495 × 39] (S3: tbl_df/tbl/data.frame)
##  $ index                      : num [1:495] 0 1 2 3 4 5 6 7 8 9 ...
##   ..- attr(*, "format.stata")= chr "%12.0g"
##  $ cz                         : num [1:495] 100 200 301 302 401 402 500 601 602 700 ...
##   ..- attr(*, "format.stata")= chr "%12.0g"
##  $ czname                     : chr [1:495] "Johnson City" "Morristown" "Middlesborough" "Knoxville" ...
##   ..- attr(*, "format.stata")= chr "%22s"
##  $ e_rank_b                   : num [1:495] 43.6 39.6 46.9 NA 43.2 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_race_bla                : num [1:495] 0.069393 0.445318 0.015122 -0.000411 -0.019341 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_race_theil_2000         : num [1:495] 0.05769 0.07497 0.11006 -0.00446 0.26349 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc               : num [1:495] 0.0409 -0.0259 0.0954 0.0105 0.0289 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc_pov25         : num [1:495] -0.0171 0.0353 0.0139 0.0545 0.0478 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc_aff75         : num [1:495] 0.02664 0.0089 0.05273 0.00763 0.06856 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_traveltime_lt15       : num [1:495] 0.384 0.354 0.479 0.661 0.406 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ hhinc00                    : num [1:495] 38210 25123 32701 36524 30403 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gini                       : num [1:495] 0.419 0.487 0.307 0.241 0.451 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ inc_share_1perc            : num [1:495] 5.7 11.82 11.81 NA 9.97 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gini99                     : num [1:495] 0.253 0.393 0.208 NA 0.348 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_middleclass           : num [1:495] 0.584 0.356 0.668 NA 0.538 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ taxrate                    : num [1:495] 0.0159 0.0257 0.0117 0.0217 0.0452 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ subcty_total_expenditure_pc: num [1:495] 1405 1569 3116 2514 2138 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ tax_st_diff_top20          : num [1:495] 1.07949 0.01987 -0.00462 0.20946 0.0028 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ eitc_exposure              : num [1:495] 0.000078 -0.012866 0.034301 1.929972 -0.003527 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ ccd_exp_tot                : num [1:495] 6.58 4.41 6.61 7.89 7.05 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ ccd_pup_tch_ratio          : num [1:495] 18 17.8 17.9 12.6 20.9 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ score_r                    : num [1:495] NA -16.02 -9.14 NA -10.45 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ dropout_r                  : num [1:495] NA 0.0562 0.0246 NA NA ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ num_inst_pc                : num [1:495] 0.0554 0.0379 0.0635 NA 0.0419 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ tuition                    : num [1:495] 1032 970 1395 NA 2524 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gradrate_r                 : num [1:495] -0.221 -0.0416 -0.1717 NA 0.1363 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_labforce                : num [1:495] 0.653 0.527 0.707 0.623 0.636 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_elf_ind_man             : num [1:495] 0.06 0.2087 0.0363 0.0685 0.1408 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ d_tradeusch_pw_1990        : num [1:495] NA 0.767 0.0275 0.0173 0.8333 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_worked1416            : num [1:495] 0.0149 0.0224 0.0275 NA 0.0239 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ mig_inflow                 : num [1:495] 0.0417 0.0182 0.031 NA 0.058 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ mig_outflow                : num [1:495] 0.0559 0.0138 0.0422 NA 0.0682 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_born_foreign            : num [1:495] 0.1605 0.0232 0.1318 0.0253 0.1903 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ scap_ski90pcm              : num [1:495] NA -0.613 -2.064 2.343 -0.694 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ rel_tot                    : num [1:495] 0.536 0.647 0.422 0.688 0.387 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ crime_violent              : num [1:495] 0.0483 0.0316 0.0512 0.0371 0.05 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_fam_wkidsinglemom       : num [1:495] 0.242 0.35 0.183 0.129 0.222 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_divorced                : num [1:495] 0.139 0.104 0.14 0.113 0.115 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_married                 : num [1:495] 0.549 0.543 0.622 0.71 0.57 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"

Le format est celui de stata, mais la plupart des variables sont numériques, cela semble correct.

# Verifier données manquantes
for (column_ in names(mobility_data)){
  print(
    paste(
      column_, 
      sum(is.na(mobility_data[, column_])), 
      sep=': '
    )
  )
}

## [1] "index: 0"
## [1] "cz: 0"
## [1] "czname: 0"
## [1] "e_rank_b: 19"
## [1] "cs_race_bla: 0"
## [1] "cs_race_theil_2000: 0"
## [1] "cs00_seg_inc: 0"
## [1] "cs00_seg_inc_pov25: 0"
## [1] "cs00_seg_inc_aff75: 0"
## [1] "frac_traveltime_lt15: 0"
## [1] "hhinc00: 0"
## [1] "gini: 0"
## [1] "inc_share_1perc: 19"
## [1] "gini99: 19"
## [1] "frac_middleclass: 19"
## [1] "taxrate: 1"
## [1] "subcty_total_expenditure_pc: 2"
## [1] "tax_st_diff_top20: 0"
## [1] "eitc_exposure: 0"
## [1] "ccd_exp_tot: 7"
## [1] "ccd_pup_tch_ratio: 20"
## [1] "score_r: 21"
## [1] "dropout_r: 89"
## [1] "num_inst_pc: 100"
## [1] "tuition: 104"
## [1] "gradrate_r: 107"
## [1] "cs_labforce: 0"
## [1] "cs_elf_ind_man: 0"
## [1] "d_tradeusch_pw_1990: 14"
## [1] "frac_worked1416: 19"
## [1] "mig_inflow: 10"
## [1] "mig_outflow: 10"
## [1] "cs_born_foreign: 0"
## [1] "scap_ski90pcm: 14"
## [1] "rel_tot: 0"
## [1] "crime_violent: 17"
## [1] "cs_fam_wkidsinglemom: 0"
## [1] "cs_divorced: 0"
## [1] "cs_married: 0"

Le nombre de points de données manquants ne semble pas non plus trop inquiétant. Il pourrait y avoir une certaine corrélation entre les points manquants mais nous pouvons ignorer cela pour le moment.

5.3 Renommer une variable

2. La variable d’intérêt est e_rank_b, renommez cette colonne. Que se passe-t-il si vous utilisez un espace ” ” dans le nom ? Est-ce une bonne idée ?

# renommer avec espace
print(names(mobility_data))

##  [1] "index"                       "cz"                         
##  [3] "czname"                      "e_rank_b"                   
##  [5] "cs_race_bla"                 "cs_race_theil_2000"         
##  [7] "cs00_seg_inc"                "cs00_seg_inc_pov25"         
##  [9] "cs00_seg_inc_aff75"          "frac_traveltime_lt15"       
## [11] "hhinc00"                     "gini"                       
## [13] "inc_share_1perc"             "gini99"                     
## [15] "frac_middleclass"            "taxrate"                    
## [17] "subcty_total_expenditure_pc" "tax_st_diff_top20"          
## [19] "eitc_exposure"               "ccd_exp_tot"                
## [21] "ccd_pup_tch_ratio"           "score_r"                    
## [23] "dropout_r"                   "num_inst_pc"                
## [25] "tuition"                     "gradrate_r"                 
## [27] "cs_labforce"                 "cs_elf_ind_man"             
## [29] "d_tradeusch_pw_1990"         "frac_worked1416"            
## [31] "mig_inflow"                  "mig_outflow"                
## [33] "cs_born_foreign"             "scap_ski90pcm"              
## [35] "rel_tot"                     "crime_violent"              
## [37] "cs_fam_wkidsinglemom"        "cs_divorced"                
## [39] "cs_married"

names(mobility_data)[4] <- "mobilite intergen"

# mobility_data$`mobilite intergen` -> On a toujours besoin des symboles `. Cela semble encombrant. 
# Essayez d'utiliser des tirets bas (_) ou des points (.)

names(mobility_data)[4] <- "mobilite_intergen"
# mobility_data$mobilite_intergen -> C'est mieux

Les réponses sont directement dans le code, moins joli mais pour cette question acceptable.

5.4 Analyse visuelle de la mobilité

3. Visualiser la distribution de la mobilité intergénérationnelle. Donnez à vos graphiques un titre approprié et des étiquettes d’axe correctes.

Nous allons analyser la densité en utilisant le lissage, un histogramme et la fonction de répartition (cdf).

hist(mobility_data$mobilite_intergen, 
     main='Intergenerational Mobility Scores', 
     xlab='Value', 
     ylab='Density', 
     probability = T, 
     breaks = seq(27, 67, 3))

lines(density(mobility_data$mobilite_intergen[!is.na(mobility_data$mobilite_intergen)]), 
      col='red', 
      lwd=2, 
      lty=2)

legend(55,0.07,
       legend=c("Density"),
       col=c("red"),
       lty=2,
       cex=0.8)

quantiles_ <- quantile(mobility_data$mobilite_intergen, 
                       probs = seq(0,1,0.01), 
                       na.rm = T)

plot(x=quantiles_, 
     y=seq(0,1,0.01), 
     type='l', 
     lwd=2, 
     main='ECDF', 
     ylab='', 
     xlab='Value for Mobility')

Les données sont principalement dispersées entre 30 et 60, il n’y a pas de sauts dans le cdf. Les données semblent également avoir un coefficient d’asymétrie positif.

5.5 Analyse visuelle: Manufacturing vs Mobilité

4. Examinez la relation entre la mobilité intergénérationnelle et la taille du secteur manufacturier. Quelle est votre conclusion ?

Si les données sont continues, commencez toujours par un nuage de points.

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

Il semble que nous ayons ici une relation linéaire. Nous pouvons également explorer cela avec une régression linéaire très simple.

# Model simple pour la relation lineaire
summary(lm(mobilite_intergen ~ cs_elf_ind_man, data=mobility_data))

## 
## Call:
## lm(formula = mobilite_intergen ~ cs_elf_ind_man, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5226  -3.6620  -0.8982   3.1683  17.9719 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     46.9024     0.5867  79.946  < 2e-16 ***
## cs_elf_ind_man -17.7064     3.1306  -5.656 2.68e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.482 on 474 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.06322,    Adjusted R-squared:  0.06125 
## F-statistic: 31.99 on 1 and 474 DF,  p-value: 2.683e-08

Les coefficients ont les signes que l’on attend. Dans ce cas, il ne semble pas que l’ajustement aux données soit excellent (R-carré très faible). Dans un cas aussi simple, c’est une bonne idée de tracer également la ligne de régression.

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

abline(a=46.9024,
       b=-17.7064,
       col='red', 
       lwd=2, 
       lty=2)

legend(0.35,60,
       legend=c("Linear Fit"),
       col=c("red"),
       lty=2,
       cex=0.8)

5.6 Analyse descriptif (secteur manufacturier)

5. Découpez la variable de la taille du secteur manufacturier en trois composantes (low-mid-high). Calculez les moyennes par groupe pour trois autres variables que vous jugez intéressantes.

Nous pouvons trouver de bons seuils en regardant le graphique ci-dessus ou en utilisant les quantiles. Voici l’approche “visuelle”. Lorsque nous coupons une variable, nous pouvons introduire des non-linéarités dans la relation (voir l’exercice suivant)

mobility_data$manuf_cut <- cut(mobility_data$cs_elf_ind_man, 
                               breaks = c(0,0.1,0.3, 1), 
                               labels = c('low', 'mid', 'high'))

mobility_data %>% 
  group_by(manuf_cut) %>% 
  summarise(
    mean_labforce = mean(cs_labforce, na.rm=T), 
    mean_growth_imports = mean(d_tradeusch_pw_1990, na.rm=T), 
    mean_mig_inflow = mean(mig_outflow, na.rm=T)
  )

## # A tibble: 3 × 4
##   manuf_cut mean_labforce mean_growth_imports mean_mig_inflow
##   <fct>             <dbl>               <dbl>           <dbl>
## 1 low               0.633               0.289          0.0458
## 2 mid               0.641               1.30           0.0427
## 3 high              0.667               2.73           0.0399

Les tendances des trois variables semblent bien correspondre aux “cuts” - mais est-ce suffisant ?

5.7 Model simple

6. Ajustez une régression qui décrit l’équation suivante \(\mathbb{E}[y| \text{manufacturing}]=\hat{\beta}_0 + \hat{\beta}_1*\text{manuf_1} + \hat{\beta}_2*\text{manuf_2}\). Pourquoi exclure la dernière composante?

R nous aide beaucoup ici:

model_lineaire_simple = lm(mobilite_intergen ~ manuf_cut, 
                           data=mobility_data)

summary(model_lineaire_simple)

## 
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut, data = mobility_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.271  -3.847  -0.658   2.952  17.605 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    45.9575     0.5388  85.293  < 2e-16 ***
## manuf_cutmid   -2.3677     0.6141  -3.856 0.000131 ***
## manuf_cuthigh  -6.3979     1.2444  -5.142    4e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.495 on 473 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.06088,    Adjusted R-squared:  0.05691 
## F-statistic: 15.33 on 2 and 473 DF,  p-value: 3.534e-07

R a automatiquement laissé tomber la première variable (pourquoi ?), le R-carré n’augmente que légèrement, ce n’était peut-être pas une si bonne idée. Nous traçons à nouveau les résultats:

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

x <- seq(0, 0.45, 0.01)
fx <- 45.9575 + (x > 0 & x <=0.1) * 0 +
      (x > 0.1 & x <= 0.3) * -2.3677 +
      (x > 0.3 & x <= 2.119) *  -6.3979
lines(x, fx, 
      col='red', 
      lty=2, 
      lwd=2)

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_simple <- predict(model_lineaire_simple, 
                              mobility_data)

points(mobility_data$cs_elf_ind_man, 
       predictions_simple, 
       col='red')

La différence entre les deux graphiques est que le plus bas montre les points réels “régressés” et le plus haut montre la ligne de régression entière.

Nous avons seulement considéré un modèle constant par morceaux, peut-être pouvons-nous obtenir de meilleurs résultats avec un modèle linéaire par morceaux ?

model_lineaire_interaction = lm(mobilite_intergen ~ manuf_cut + cs_elf_ind_man + manuf_cut*cs_elf_ind_man, 
                           data=mobility_data)

summary(model_lineaire_interaction)

## 
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut + cs_elf_ind_man + 
##     manuf_cut * cs_elf_ind_man, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.9773  -3.6249  -0.6577   2.7890  17.3450 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    42.921      1.983  21.640   <2e-16 ***
## manuf_cutmid                    3.175      2.252   1.410   0.1592    
## manuf_cuthigh                  -5.908     10.728  -0.551   0.5821    
## cs_elf_ind_man                 44.485     27.973   1.590   0.1124    
## manuf_cutmid:cs_elf_ind_man   -57.912     28.507  -2.031   0.0428 *  
## manuf_cuthigh:cs_elf_ind_man  -37.336     40.598  -0.920   0.3582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.463 on 470 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.07768,    Adjusted R-squared:  0.06787 
## F-statistic: 7.917 on 5 and 470 DF,  p-value: 3.562e-07

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_interact <- predict(model_lineaire_interaction, 
                              mobility_data)

points(mobility_data$cs_elf_ind_man, 
       predictions_interact, 
       col='red')

Nous pouvons constater que la fraction de travailleurs manufacturiers n’a pas un effet régulier sur la mobilité intergénérationnelle. Cela pose généralement un problème lorsque nous travaillons avec des modèles linéaires. Nous avons vu ici une possibilité de résoudre le problème. En subdivisant la variable explicative, nous pouvons obtenir des non-linéarités (nous y reviendrons plus tard dans le cours) dans un modele lineaire.

5.8 Model “Classe moyenne”

L’existence d’une classe moyenne importante est souvent considérée comme une bonne chose. Voir par exemple cet article, analysons cela plus en détail.

7. Réalisez un nuage de points (“scatterplot”) avec la mobilité intergénérationnelle sur l’axe des y et la fraction de la classe moyenne sur l’axe des x. Quelles sont vos conclusions ?

8. Exécutez un modèle linéaire qui explique la mobilité intergénérationnelle à l’aide de la fraction de la classe moyenne.

plot(x=mobility_data$frac_middleclass, 
     y=mobility_data$mobilite_intergen, 
     main='Middle Class and Intergen. Mobility', 
     ylab='Intergen. Mobility', 
     xlab='Fraction Middle class')

Déjà dans le graphique, nous pouvons voir une relation entre les deux variables, continuons avec une régression.

model_middle_class <- lm(mobilite_intergen ~ frac_middleclass, 
                         data=mobility_data)

summary(model_middle_class)

## 
## Call:
## lm(formula = mobilite_intergen ~ frac_middleclass, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.9940  -2.9766  -0.5619   2.7327  15.9744 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        18.986      1.381   13.74   <2e-16 ***
## frac_middleclass   45.411      2.491   18.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.343 on 474 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.4121, Adjusted R-squared:  0.4109 
## F-statistic: 332.3 on 1 and 474 DF,  p-value: < 2.2e-16

plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_middleclass <- predict(model_middle_class, 
                              mobility_data)

points(mobility_data$frac_middleclass,
       predictions_middleclass, 
       col='red')

La régression semble être en accord avec l’article, mais attention! Nous n’avons considéré que la corrélation entre les deux variables (correlation \(\neq\) causation). Il y a beaucoup de possibilités qui pourraient générer les données. Par exemple, il pourrait y avoir une troisième variable qui a un effect causal sur les deux variables que nous considérons ici.

Mais laissons cela de côté pour le moment : Après chaque régression, nous devrions toujours considérer les résidus:

# Analyse des residues
residuals_ <- mobility_data$mobilite_intergen - predictions_middleclass

# KDE smoothing (à voir plus tard)
k1_ <- mobility_data$frac_middleclass[!is.na(mobility_data$frac_middleclass)]
k2_ <- residuals_[!is.na(mobility_data$frac_middleclass)]

smooth_fit <- lowess(k1_, k2_)

plot(mobility_data$frac_middleclass, 
     residuals_, 
     ylab='Residuals', 
     xlab='Middle class')

lines(smooth_fit, lty=2, lwd=2, col='red')

Il est clair qu’il reste des informations dans les données. En utilisant le graphique ci-dessus, nous pouvons décider d’ajouter une fonction polynomiale dans les données (une autre façon pour générer les non-linéarités).

mobility_data_nomissings <- mobility_data[!is.na(mobility_data$frac_middleclass), ]

model_middle_class_poly <- lm(mobilite_intergen ~ poly(frac_middleclass, 2), 
                              data=mobility_data_nomissings)

summary(model_middle_class_poly)

## 
## Call:
## lm(formula = mobilite_intergen ~ poly(frac_middleclass, 2), data = mobility_data_nomissings)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0057  -2.9207  -0.3494   2.4430  15.5096 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 43.9039     0.1881 233.464  < 2e-16 ***
## poly(frac_middleclass, 2)1  79.1684     4.1029  19.296  < 2e-16 ***
## poly(frac_middleclass, 2)2  31.2715     4.1029   7.622 1.38e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.103 on 473 degrees of freedom
## Multiple R-squared:  0.4764, Adjusted R-squared:  0.4742 
## F-statistic: 215.2 on 2 and 473 DF,  p-value: < 2.2e-16

plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_middleclass_poly <- predict(model_middle_class_poly,
                                        mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])

lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
       predictions_middleclass_poly, 
       col='red')

L’ajustement est clairement meilleur, ici l’ajout du terme polynomial semble ajouter beaucoup de valeur (comparez eg. les R-carrés). Mais que se passe-t-il si nous ajoutons de plus en plus de termes polynomiaux ?

5.9 Challenge

# List pour les resultats
results_list <- list()

# Loop
for (length_ in seq(1,10,1)){
  reg_ <- lm(mobilite_intergen ~ poly(frac_middleclass, length_), 
             data=mobility_data_nomissings)
  
  print(paste('R-Squared: ', summary(reg_)$r.squared))
  
  predictions_ <- predict(reg_,
                          mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])
  
  results_list[[length_]] <- predictions_

}

## [1] "R-Squared:  0.412134418340672"
## [1] "R-Squared:  0.476437415437812"
## [1] "R-Squared:  0.477700187650047"
## [1] "R-Squared:  0.477727823275877"
## [1] "R-Squared:  0.478923076244006"
## [1] "R-Squared:  0.48017980394324"
## [1] "R-Squared:  0.48426795442976"
## [1] "R-Squared:  0.487546350228303"
## [1] "R-Squared:  0.487606950061351"
## [1] "R-Squared:  0.489193655479349"

plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

for(i in seq(1,10,1)){
  
  predictions_middleclass_poly <- results_list[[i]]
  
  lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
         predictions_middleclass_poly, 
         col='red')
  
}

Le R-carré augmente toujours avec des variables supplémentaires, mais cela ne signifie pas nécessairement que le modèle est meilleur. Vous devez toujours vérifier le modèle et l’ajuster en conséquence.